/*---------------------------------------------------------------------------*\
License
    This file is not part of OpenFOAM but a modified version of original 
    routine distributed with OpenFOAM.

    You can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is a trademark of OpenCFD.

    See GNU General Public License at <http://www.gnu.org/licenses/>

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    rough wall function for mut

Description
    ATS4i fixed original code bugs regarding SIGFPE
    All equations are the same but with correct implementation

    This version is not endorsed by original OpenFoam authors. Use at your own risk.

Authors
    Original code by OpenCFD
    Modified version by ATS4i - Dec/2010

\*---------------------------------------------------------------------------*/

#include "mutRoughWallFunctionFvPatchScalarField.H"
#include "RASModel.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace compressible
{
namespace RASModels
{


// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

scalar mutRoughWallFunctionFvPatchScalarField::fnRough
(
    const scalar KsPlus,
    const scalar Cs
) const
{
    // Return fn based on non-dimensional roughness height

    if (KsPlus < 90.0)
    {
        return pow
        (
            (KsPlus - 3.328)/86.672 + Cs*KsPlus,
            sin(0.4258*(log(max(1e-14,KsPlus)) - 0.811))
        );
    }
    else
    {
        return (1.0 + Cs*KsPlus);
    }
}


tmp<scalarField> mutRoughWallFunctionFvPatchScalarField::calcMut() const
{
    const label patchI = patch().index();

    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
    const scalarField& y = rasModel.y()[patchI];
    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
    const tmp<volScalarField> tk = rasModel.k();
    const volScalarField& k = tk();
    const scalarField& muw = rasModel.mu().boundaryField()[patchI];

    const scalar Cmu25 = pow(Cmu_, 0.25);

    tmp<scalarField> tmutw(new scalarField(*this));
    scalarField& mutw = tmutw();

    forAll(mutw, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        scalar uStar = Cmu25*sqrt(k[faceCellI]);

        scalar yPlus = uStar*y[faceI]/(muw[faceI]/rhow[faceI]);

        scalar KsPlus = uStar*Ks_[faceI]/(muw[faceI]/rhow[faceI]);

	scalar Edash = E_;
//        scalar Edash = 0;
        scalar yPlusLamNew = yPlusLam_;
//	scalar CsKays = max(1e-14, 0.3005213-1/KsPlus);

        if (KsPlus > 3.328) //(32.6/9.8)) //2.25)
        {
	    /* it was
            yPlusLamNew = rasModel.yPlusLam(kappa_, max(Edash,1e-14));
            Edash /= fnRough(KsPlus, Cs_[faceI]); */
	    
	    /* it becomes */
            Edash = 1/kappa_*log(fnRough(KsPlus, Cs_[faceI]));

	    //Edash = 1/kappa_*log(fnRough(KsPlus, CsKays));
            yPlusLamNew = rasModel.yPlusLam(kappa_, max(E_,1e-14));
        }

        /*if (debug)
        {
            Info<< "yPlus = " << yPlus
                << ", KsPlus = " << KsPlus
                << ", Edash = " << Edash
                << ", yPlusLam = " << yPlusLam_
                << endl;
        }*/

        if (yPlus > yPlusLamNew)
        {
	    /* it was
            mutw[faceI] =
                muw[faceI]*(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1);

	    Calculation from OpenFOAM original mutWallFunction (smooth wall!!)
            mutw[faceI] = muw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1);*/

            mutw[faceI] = muw[faceI]*(yPlus*kappa_/
	    log(max(E_*yPlus/exp(kappa_*Edash),1+1e-4)) - 1);

        }
    }

    return tmutw;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    mutWallFunctionFvPatchScalarField(p, iF),
    Ks_(p.size(), 0.0),
    Cs_(p.size(), 0.0)
{}


mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
    const mutRoughWallFunctionFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    mutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
    Ks_(ptf.Ks_, mapper),
    Cs_(ptf.Cs_, mapper)
{}


mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    mutWallFunctionFvPatchScalarField(p, iF, dict),
    Ks_("Ks", dict, p.size()),
    Cs_("Cs", dict, p.size())
{}


mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
    const mutRoughWallFunctionFvPatchScalarField& rwfpsf
)
:
    mutWallFunctionFvPatchScalarField(rwfpsf),
    Ks_(rwfpsf.Ks_),
    Cs_(rwfpsf.Cs_)
{}


mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
    const mutRoughWallFunctionFvPatchScalarField& rwfpsf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    mutWallFunctionFvPatchScalarField(rwfpsf, iF),
    Ks_(rwfpsf.Ks_),
    Cs_(rwfpsf.Cs_)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void mutRoughWallFunctionFvPatchScalarField::autoMap
(
    const fvPatchFieldMapper& m
)
{
    mutWallFunctionFvPatchScalarField::autoMap(m);
    Ks_.autoMap(m);
    Cs_.autoMap(m);
}


void mutRoughWallFunctionFvPatchScalarField::rmap
(
    const fvPatchScalarField& ptf,
    const labelList& addr
)
{
    mutWallFunctionFvPatchScalarField::rmap(ptf, addr);

    const mutRoughWallFunctionFvPatchScalarField& nrwfpsf =
        refCast<const mutRoughWallFunctionFvPatchScalarField>(ptf);

    Cs_.rmap(nrwfpsf.Cs_, addr);
    Ks_.rmap(nrwfpsf.Ks_, addr);
}


void mutRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
{
    fvPatchField<scalar>::write(os);
    writeLocalEntries(os);
    Cs_.writeEntry("Cs", os);
    Ks_.writeEntry("Ks", os);
    writeEntry("value", os);
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

makePatchTypeField(fvPatchScalarField, mutRoughWallFunctionFvPatchScalarField);

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace RASModels
} // End namespace compressible
} // End namespace Foam

// ************************************************************************* //
